c
c-----------------------------------------------------------------------------
c
c  Naca6.f -- contains routines to generate the NACA 1 and 6 series
c             foils.
c
c  Written By: S.E.Norris
c
c  norris@cfd.mech.unsw.edu.au
c
c  $RCSfile: naca6.f,v $
c  $Author: norris $
c  $Revision: 1.5 $
c  $Date: 1995/08/31 11:05:52 $
c
c  $Log: naca6.f,v $
c  Revision 1.5  1995/08/31  11:05:52  norris
c  *** empty log message ***
c
c  Revision 1.4  1995/08/31  07:30:19  norris
c  Cleaned up code.
c
c  Revision 1.3  1995/08/31  03:35:01  norris
c  Changed to use the directory.h header file that holds the directory name.
c
c  Revision 1.2  1995/08/31  02:03:05  norris
c  Added the use of the lnblnk function to generate pathnames.
c
c  Revision 1.1  1995/08/30  11:07:20  norris
c  Initial revision
c
c-----------------------------------------------------------------------------
c
      SUBROUTINE Naca1( x,npl,npmx,naca,inaca,scle,t_n )
c
      IMPLICIT none
      INTEGER npl,npmx,inaca
      REAL x(3,npmx),scle
      CHARACTER naca*(*)
      LOGICAL t_n
c
c  This routine calculates the ordinates of a NACA 1 series foil.
c
      INTEGER nbmx
      PARAMETER (nbmx=30)
      REAL pi
      PARAMETER (pi=3.14159265358979)
c
      INTEGER npt,npf,nb,nbm
      INTEGER i,j
      REAL ndata(3)
      REAL ab,bb,rb,xb(nbmx),yb(nbmx),y2b(nbmx)
      REAL xi,xm
      REAL dx,dy,dydx,t,y
      CHARACTER filnme*6
      REAL Nacat
      EXTERNAL Nacat
c
      DATA  filnme / '16-015' /
c
      npt  = npl+1
      npf  = npl/2
c
c     Calculate the parameters describing the foil, from the given NACA
c     number. Then calculate the cubic spline for the camber.
c
      call Naca1b( naca,inaca,ndata )
      call GetSpl( xb,yb,y2b,nb,nbm,nbmx,ab,bb,rb,filnme,t_n,ndata(3) )
      if (.not. t_n) return
c
c     Set the x values according to a cosine distribution.
c
      call Cosdn( x,npl,npmx )
c
c     Calculate 'y' values for each value of 'x'.
c
      x(2,1)     = 0.0
      x(2,npt)   = 0.0
      x(2,npf+1) = 0.0
      do i = 2,npf
	 j = npt+1-i
c
c        Calculate the camber line.
c
	 xi = x(1,i)
	 xm = 1.0 - x(1,i)
	 if (xi.le.0.0 .or. xi.ge.1.0) then
	    y    = 0.0
	    dydx = 0.0
	 else
	    y    = (-ndata(1)/( 4.0*pi))*( xm*log(xm) + xi*log(xi) )
	    dydx = ( ndata(1)/( 4.0*pi))*(    log(xm) -    log(xi) )
	 endif
c
c        Calculate thickness.
c
	 t = Nacat( xi,xb,yb,y2b,nbm,nb,ab,bb )
c
c        Finally, calculate node positions upon the foil surface.
c
	 dy = sqrt(( t**2 )/( 1.0+( dydx**2 )))
	 dx = dy * dydx
c
	 x(1,i) = x(1,i) + dx
	 x(1,j) = x(1,i) - 2.0*dx
	 x(2,i) = y - dy
	 x(2,j) = y + dy
      enddo
c
c     Rescale dimensions if nessisary.
c
      call FoilScale( x,npt,scle )
c
      return
      END
c
c-----------------------------------------------------------------------------
c
      SUBROUTINE Naca1b( naca,inaca,ndata )
c
      IMPLICIT none
      INTEGER inaca
      REAL ndata(3)
      CHARACTER naca*(*)
c
c  This routine converts a NACA number into the parameters describing
c  a foil.
c
      INTEGER icld,itau
      INTEGER Ist
      EXTERNAL Ist
c
c     Convert first four digits into integers.
c
      inaca= Ist( 2,naca(5:6),0 )
      icld = Ist( 1,naca(4:4),0 )
      itau = Ist( 2,naca(5:6),0 )
      ndata(1) =   0.1*real(icld)
      ndata(2) =  0.01*real(itau)
      ndata(3) =  real(itau)/15.0
c
      return
      END
c
c-----------------------------------------------------------------------------
c
      SUBROUTINE Naca6( x,npl,npmx,naca,inaca,scle,a,type,t_n )
c
      IMPLICIT none
      INTEGER npl,npmx,inaca
      REAL x(3,npmx),scle,a
      CHARACTER naca*(*),type*3
      LOGICAL t_n
c
c  This routine calculates the ordinates of a NACA 6 series foil.
c
      INTEGER nbmx
      PARAMETER (nbmx=30)
c
      INTEGER npt,npf,nb,nbm
      INTEGER i,j
      REAL ndata(3),dx,dy,dydx,t,y
      REAL ab,bb,rb,xb(nbmx),yb(nbmx),y2b(nbmx)
      CHARACTER filnme*6
      REAL Nacat
      EXTERNAL Nacat
c
      npt  = npl+1
      npf  = npl/2
c
c     Calculate the parameters describing the foil, from the given NACA
c     number. Then calculate the cubic spline for the camber.
c
      call Naca6b( naca,inaca,ndata,filnme,type )
      if (ndata(2).ne.0.0) then
	 call GetSpl(xb,yb,y2b,nb,nbm,nbmx,ab,bb,rb,filnme,t_n,ndata(3))
      endif
      if (.not. t_n) return
c
c     Set the x values according to a cosine distribution.
c
      call Cosdn( x,npl,npmx )
c
c     Calculate 'y' values for each value of 'x'.
c
      x(2,1)     = 0.0
      x(2,npt)   = 0.0
      x(2,npf+1) = 0.0
      do i = 2,npf
	 j = npt+1-i
c
c        Calculate the camber line.
c
	 call Naca6c( y,dydx,x(1,i),a,ndata(1) )
c
c        Calculate thickness.
c
	 if (ndata(2).eq.0.0) then
	    t = 0.0
	 else
	    t = Nacat( x(1,i),xb,yb,y2b,nbm,nb,ab,bb )
	 endif
c
c        Finally, calculate node positions upon the foil surface.
c
	 dy = sqrt(( t**2 )/( 1.0+( dydx**2 )))
	 dx = dy * dydx
c
	 x(1,i) = x(1,i) + dx
	 x(1,j) = x(1,i) - 2.0*dx
	 x(2,i) = y - dy
	 x(2,j) = y + dy
      enddo
c
c     Rescale dimensions if nessisary.
c
      call FoilScale( x,npt,scle )
c
      return
      END
c
c-----------------------------------------------------------------------------
c
      SUBROUTINE Naca6b( naca,inaca,ndata,filnme,type )
c
      IMPLICIT none
      INTEGER inaca
      REAL ndata(3)
      CHARACTER naca*(*),filnme*6,type*3
c
c  This routine converts a NACA number into the parameters describing
c  a foil.
c
      INTEGER icld,itau,ibase,i
      INTEGER Ist
      EXTERNAL Ist
c
c
      if (type(2:2).ne.'(') then
	 filnme = naca(1:3)//'0'//naca(5:6)
	 ibase  = Ist( 2,naca(5:6),0 )
	 i = -3
      else
	 if (type(3:3).eq.'1') then
	    i = 0
	    filnme = naca(1:2)//type(1:1)//'00'//naca(4:4)
	 else if (type(3:3).eq.'2') then
	    i = 1
	    filnme = naca(1:2)//type(1:1)//'0'//naca(4:5)
	 endif
	 ibase  = Ist( 2,naca(4:4+i),0 )
      endif
c
      inaca= Ist( 2,naca(8+i:9+i),0 )
      icld = Ist( 1,naca(7+i:7+i),0 )
      itau = Ist( 2,naca(8+i:9+i),0 )
      ndata(1) =  0.1*real(icld)
      ndata(2) = 0.01*real(itau)
      ndata(3) = real(itau)/real(ibase)
c
      return
      END
c
c-----------------------------------------------------------------------------
c
      SUBROUTINE Naca6c( y,dydx,x,a,cl )
c
      IMPLICIT none
      REAL y,dydx,x,a,cl
c
c  Calculates the camber line for 6-series foils
c
      REAL pi
      PARAMETER (pi=3.14159265358979)
      REAL xm,am,g,h
      REAL ax,da,ya
c
      xm = 1.0-x
      am = 1.0-a
      ax = a-x
      if (x.le.0.0 .or. x.ge.1.0) then
	 y    = 0.0
	 dydx = 0.0
      else
	 if (a.eq.1.0) then
	    y    = (-cl/( 4.0*pi))*( xm*log(xm) + x*log(x) )
	    dydx = ( cl/( 4.0*pi))*(    log(xm) -   log(x) )
	 else
	    if (a.eq.0.0) then
	       g = -0.25
	    else
	       g = (-1.0/am )*( a*a*( 0.5*log(a ) - 0.25 ) + 0.25 )
	    endif
	    h = ( 1.0/am )*( 0.5*am*am*log(am) - 0.25*am*am ) + g
	    if (x.ne.a) then
	       ya = ( 0.5*ax*ax*log(abs(ax)) - 0.5*xm*xm*log(xm)
     &            +  0.25*xm*xm - 0.25*ax*ax )/am - x*log(x) + g - h*x
	       da = (    log(xm)*( xm + 0.5*xm*xm ) + 0.5*xm*xm
     &            - log(abs(ax))*( ax + 0.5*ax*ax ) - 0.5*ax*ax
     &            - 0.5*xm + 0.5*ax )/am - log(x) - h - 1.0
	    else
	       ya = ( -0.5*xm*xm*log(xm) + 0.25*xm*xm )/am
     &            - x*log(x) + g - h*x
	       da = (    log(xm)*( xm + 0.5*xm*xm ) + 0.5*xm*xm
     &            - 0.5*xm )/am - log(x) - h - 1.0
	    endif
	    y    = cl*ya/( 2.0*pi*( 1.0+a ) )
	    dydx = cl*da/( 2.0*pi*( 1.0+a ) )
	 endif
      endif
c
      return
      END
c
c-----------------------------------------------------------------------------
c
      SUBROUTINE GetSpl( xb,yb,y2b,nb,nbm,nbmx,ab,bb,rb,filnme,t_n,scl )
c
      IMPLICIT none
      INTEGER nb,nbm,nbmx
      REAL xb(nbmx),yb(nbmx),y2b(nbmx),ab,bb,rb,scl
      CHARACTER filnme*(*)
      LOGICAL t_n
c
c  Calculates the thickness spline for 1-series and 6-series foils
c
      INTEGER i
      REAL y1,y2,xb2,xb3
c
      DATA  y2   / 3.0e35 /
c
c
      call GetCor( xb,yb,nb,nbmx,rb,filnme,t_n )
      if (t_n) then
	 nbm = nb-1
	 if (scl.ne.1.0) then
	    do i = 1,nb
	       yb(i) = yb(i)*scl
	    enddo
	 endif
	 xb2 = sqrt( xb(2) )
	 xb3 = sqrt( xb(3) )
	 ab  = ( yb(2)*xb(3)-yb(3)*xb(2) )/( xb(3)*xb2 - xb(2)*xb3 )
	 bb  = ( yb(2) - ab*xb2 )/xb(2)
	 y1  = bb + 0.5*ab/xb2
	 call Spln( xb(2),yb(2),nbm,y1,y2,y2b(2) )
      endif
c
      return
      END
c
c-----------------------------------------------------------------------------
c
      SUBROUTINE GetCor( xb,yb,nb,nbmx,rb,filnme,t_n )
c
      IMPLICIT none
      INTEGER nbmx,nb
      REAL xb(nbmx),yb(nbmx),rb
      CHARACTER filnme*6
      LOGICAL t_n
c
c  Reads in the coordinates of a NACA thickness distribution from file
c
      INCLUDE 'directory.h'
      INTEGER ierr,i
      REAL x,y
      CHARACTER fullpath*128
      LOGICAL fex
      INTEGER Lnblnk
      EXTERNAl Lnblnk
c
      fullpath = direct//filnme
      inquire(file=fullpath,exist=fex)
      if (.not.fex) then
         print '(2a)', 'GetCor(): Cannot find file ',fullpath
         t_n = .false.
         return
      endif
      open( unit=10,file=fullpath,status='OLD',iostat=ierr )
      if (ierr.ne.0) then
         print '(2a)', 'GetCor(): Error opening file ',fullpath
         t_n = .false.
         return
      endif
c
      read(10,*)
      read(10,*)
      x = 0.0
      i = 0
      do while( x.lt.100.0 )
	 read(10,*) x,y
	 i = i+1
	 xb(i) = 0.01*x
	 yb(i) = 0.01*y
      enddo
      read(10,*) rb
      rb = 0.01*rb
      nb = i
c
      close(unit=10)
c
      return
      END
c
c-----------------------------------------------------------------------------
c
      FUNCTION Nacat( x,xb,yb,y2b,nbm,nb,ab,bb )
c
      IMPLICIT none
      INTEGER nb,nbm
      REAL x,xb(nb),yb(nb),y2b(nb),ab,bb
      REAL Nacat
c
c  Returns an interpolated thickness for naca 1- and 6-series foils
c
      REAL y
c
      if (x.lt.xb(2)) then
	 y = ab*sqrt(x) + bb*x
      else
	 call Spin( xb(2),yb(2),y2b(2),nbm,x,y )
      endif
      Nacat = y
c
      return
      END
